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1.0  INTRODUCTION 

The  objective  of  this  research  effort  was  to  further 
investigate  the  applicability  of  statistical  signal  processing 
techniques  developed  during  the  previous  contract  (Contract  # 
F0860b-74-C-0020)  for  the  purpose  of  enhancing  seismic  source 
depth  determinations.  These  techniques  were  designed  to  improve 
the  detection  and  interpretation  of  seismic  depth  phase  information 
contained  throughout  the  seismogram.  Three  seismic  events  were 
analyzed  during  this  contract  period,  and  the  results  obtained 
clearly  demonstrated  the  ability  of  these  new  depth  phase  detec- 
tion techniques  to  obtain  seismic  source  depths  for  events  in 
which  conventional  analysis  was  ineffective. 

The  primary  interest  in  improving  seismic  depth  phase  detec- 
tion is  that  depth  phase  information  gives  reliable  estimates 
of  seismic  source  depths.  Source  depth  is  a critical  factor  in 
the  discrimination  between  man-made  events  and  natural  seismic 
events  since  the  maximum  depth  of  burial  for  man-made  events 
is  of  the  order  of  ten  kilometers  and  earthquakes  can  originate 
at  much  greater  depths . 

At  present,  depth  phase  detection  is  performed  by  analyzing 
the  first  arrival  portion  of  the  seismogram  for  the  presence 
of  the  pP  and  sP  arrivals.  Analysis  is  usually  done  through 
attempts  to  visually  identify  the  pP  and  sP  first  movements 
from  seismograms  recorded  over  a suite  of  stations.  Identifica- 
tion of  the  depth  phase  is  verified  when  the  variations  in  the 
differential  travel  times  are  in  agreement  with  the  differential 
travel  time  tables  for  the  given  source  to  station  differences. 
Correlation  techniques  are  also  applied  to  first  arrival 
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portions  jf  seismograms  in  an  effort  to  detect  the  presence 
of  the  depth  phases  and  their  delay  times.  Both  of  these 
detection  methods  use  only  the  first  arrival  portion  of  the 
seismogram  containing  the  pP  and  sP  first  movements  and  ignore 
the  abundance  of  depth  phase  information  contained  in  the  direct 
and  surface  reflected  waveforms  present  throughout  the  entire 
seismogram.  Ignoring  this  additional  depth  information  limits 
the  inherent  detectability  of  these  methods. 

Our  efforts  have  focused  on  the  development  of  a cepstral 
estimation  procedure  which  makes  constructive  use  of  the  depth 
phase  information  present  throughout  the  seismograms  recorded 
over  a suite  of  stations.  This  analysis  is  based  on  techniques 
which  account  for  the  depth  phase  delay  time  variations  asso- 
ciated with  the  later  seismic  arrivals  and  on  a technique  which 
utilizes  the  phase  differences  between  the  direct  and  surface 
reflected  arrivals  in  order  to  enhance  the  detection  of  depth 
phases.  A Cepstrum  Matched  Filter  technique  was  developed  to 
aid  in  interpreting  the  typically  complicated  cepstral  patterns 
resulting  from  this  analysis.  It  has  proven  to  be  very  valuable 
ii  further  improving  seismic  depth  phase  detection. 

Based  on  the  work  done  to  date,  all  indications  are  that 
these  new  seismic  depth  phase  detection  techniques  will  greatly 
enhance  the  analyst's  ability  to  determine  seismic  source  depths. 
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2.0  DEPTH  PHASE  ANALYSIS  TECHNIQUES 


The  most  accurate  estimates  of  teleseismic  source  depths 
are  obtained  through  the  measurement  of  the  differential  travel 
times  between  the  direct  seismic  arrival  (P  wave)  and  the 
surface  reflected  arrivals  (pP  and  sP) . The  pP  phase  follows  the 
P phase  by  roughly  1/4  second  for  each  kilometer  of  source 
depth,  and  the  sP  phase  follows  the  P phase  by  about  1/3  second 
per  kilometer  of  depth.  Therefore,  by  obtaining  estimates 
of  P-pP  or  P-sP  delay  times  to  within  1/3  second,  it  is 
possible  to  estimate  the  source  depth  to  within  one  or  two 
kilometers.  Depth  estimates  of  this  accuracy  are  extremely 
useful  in  distinguishing  natural  events  from  man-made  events. 

The  present  difficulty  in  using  depth  phase  information  to  deter- 
mine source  depths  is  that  the  depth  phases  are  identified  for 

only  about  half  the  events  analyzed  using  the  conventional  tech- 
niques . 


Depth  phase  detection  is  generally  performed  by  visual 
identification  of  the  P,  pP,  and  sP  first  movements  (See  Figure 
1) . Identifications  of  the  depth  phases  are  verified  when  the 
variations  in  the  differential  travel  times  are  in  agreement 
with  the  differential  travel  time  tables  for  the  given  source 
to  station  distances.  Both  the  detectability  and  accuracy  of 
these  delay  time  estimates  are  limited  in  this  procedure  since 
both  are  based  on  identifying  the  P,  pP,  or  sP  first  movements. 
Correlation  techniques  are  sometimes  used  to  analyze  the  first 
arrival  portion  of  the  seismogram  (that  portion  containing  the 
P,  pP,  and  sP  first  arrivals)  for  the  presence  of  the  pP  and 
sP  phases.  However,  both  of  these  detection  techniques  have 
limited  detectability  and  resolution  since  they  use  depth 
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Figure  1.  An  Example  of  Unfiltercd  Z-Component 
Seismograms  showing  Initial  Phases 
(Bolivia  Event  12:02:22  Z 19  July  1962). 
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phase  information  contained  only  in  the  first  arrival  portion 
of  the  seismogram  and  ignore  the  abundance  of  depth  phase 
information  contained  in  the  direct  and  surface  reflected 
waveforms  present  throughout  the  rest  of  the  seismogram. 

In  the  following  section,  we  discuss  statistical  techniques 
for  utilizing  depth  phase  information  contained  in  the  entire 
seismogram.  These  techniques  are  Stochastic  Cepstrum  Stacking 
and  Stochastic  Cesptrum  Phaser  Stacking. 


f 
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2 . 1 Statistical  Techniques  for  Utilizing  Depth  Phase 
Information  Contained  in  the  Coda 


In  order  to  constructively  use  the  depth  phase  information 
contained  in  the  seismic  coda,  one  must  be  aware  of  variations 
in  the  depth  phase  delay  times  observed  along  the  coda.  The 
largest  of  these  variations  originate  from  the  differential 
travel  times  between  later  seismic  arrivals  such  as  PP,  PcP, 

PPP,  etc  and  their  associated  surface  reflections  pPP,  sPP, 
pPcP,  sPcP,  pPPP,  sPFP,  etc.  Other  contributions  to  the  delay 
time  variations  include  background  noise,  multipath  effects,  and 
differential  distortions  in  the  seismic  waveforms  caused  by  the 
properties  of  the  medium.  A cepstrum  estimation  procedure 
that  allows  for  these  delay  time  variations  has  proven  to 
be  very  useful.  As  described  in  Appendix  A,  "Stochastic  Cepstrum 
Stacking,"  this  procedure  allows  for  these  delay  time  variations 
by  stochastically  averaging  cepstra  calculated  from  different 
portions  of  the  seismogram  to  enable  improved  detection  of  pP-P 
time  delays.  With  conventidnal  averaging  of  these  cepstra,  detec- 
tability is  significantly  reduced. 

In  many  cases,  further  enhancement  in  depth  phase  detecta- 
bility is  achieved  by  utilizing  the  relative  phase  difference 
between  the  direct  and  surface  reflected  seismic  arrivals.  In 
Appendix  B,  "Stochastic  Cepstrum  Phasor  Stack,"  we  show  that 
the  phase  of  the  cosinusoidal  modulation  of  the  power  spectrum 
resulting  from  the  presence  of  the  delayed  surface  reflected 
arrival  can  be  assigned  to  each  cepstrum  amplitude  and  these 
values  can  then  be  treated  as  cepstrum  phasors.  These  phase 
angles  can  be  used  to  help  distinguish  cepstrum  peaks  resulting 
from  surface  reflected  arrivals  from  those  resulting  from  spu- 
rious sources. 
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Both  the  Stochastic  Cepstrum  Stack  and  Stochastic  Cepstrum 
Phasor  Stack  techniques  were  used  in  the  analysis  of  events 
during  this  contract.  The  results  are  given  in  Section  3.  The 
detailed  mathematical  description  of  the  steps  involved  in  com- 
puting these  cepstra  is  described  in  Appendix  C. 

In  the  next  section,  we  discuss  a Cepstrum  Matched  Filter 
technique  which  was  developed  as  an  aid  to  interpretating  the 
depth  phc'se  information  contained  in  these  computed  cepstra. 


t 
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2-2  Cepstrum  Matched  Filter  Technique 


\ 


The  previous  sections  described  techniques  which  enable 
one  to  utilize  the  depth  phase  information  contained  in  the 
seism. c coda  and  thereby  enhance  the  appearance  of  the  cepstrum 
peak  corresponding  to  the  pP-P  delay  time.  If  one  can  identify 
such  a cepstrum  peak,  the  source  depth  can  then  be  accurately  ’ 
determined  from  differential  travel  time  tables,  provided  the 
velocity  versus  depth  assumptions  are  valid.  However,  cepstrum 
peaks  corresponding  to  the  pP-P  time  delay  Tj,  the  sP-P  time 
delay  r2,  as  well  as  the  sP-pP  time  delay,  will  In  general  result 
from  a single  seismic  source.  According  th  the  particular  cepstru, 
algorithm  used,  additional  peaks  will  appear  at  nt  ♦ mr  (where 
n and  m arc  integers)  when  non-linear  operations  are  included  in 
the  cepstrum  computation.  Thus,  there  is  in  general  a complicated 
set  of  cepstral  peaks  corresponding  to  a single  source  depth. 

In  this  section,  we  describe  the  Cepstrum  Matched  Filter  (CMF) 
technique,  which  is  an  automatic  method  of  recogniting  the  presence 
of  a cepstrum  pattern  corresponding  to  a particular  source  depth, 
over  a range  of  possible  source  depths. 

If  one  assumes  that  the  cepstral  peaks  appearing  at  the 

PP"P’  SP'P  and  SP~pP  time  delays  dominate  the  expected  cepstrum 
pattern  for  an  event,  then  one  can  formulate  the  CMF  algorithm 
m the  following  way.  Consider  the  cepstrum  to  consist  of  N 
amplitude  values  CP(n)  for  the  time  delays  tn= (n- 1) • At , where 
n = 1 , 2 , 3 . . . N . Then  the  synthetic  cepstrum  (CS)  expected  for  an 
event  at  a given  depth  and  distance  corresponding  to  a pP-P  delay 
time  of  t,  can  be  represented  by: 
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CS (n , i , a)  = Q (n , (a  - 1) t ) + Q(n,x)  + Q(n,ax) 


with  a-  Tsp.p/Tpp.p»  and  is  the  ratio  of  delay  times  for  sP-P 
and  pP-P.  Also  Q is  defined  by 

Q(n,x)  = 1 ; if  x-At<tn<x+At  , 

= 0 ; otherwise  . 


Thus,  CS  is  represented  by  three  unit  amplitude  peaks  located  at 
the  sP-pP,  pP-P,  and  sP-P  delay  times,  each  peak  consisting  of 
three  adjacent  delay  time  points. 


The  CMP  output  at  time  delay  t is  defined  to  be  the  maximum 
zero  lag  correlation  of  the  computed  ccpstrum  (CP)  with  the 
synthesized  ccpstrum  (CS)  for  a source  d’pth  having  a pP-P  delay 
time  t computed  over  a range  of  ratios  «.  This  can  be  written 
as 


CMF(t) 


Max 

otj_<  a <_  ou, 


N 

I 

n=  1 


CP (n) *CS(n,i 


if  CP (x ) > . 7 *CP  (ax ) 
and  CP(x)>.7*CP((a-l) r) 


= CP(x)  ; otherwise 


Here,  aj  and  02  are  set  from  the  expected  range  of  values 
for  a for  a given  range  of  possible  depths.  Values  of  Cj  = 1.25 
and  a2  = 1.55  were  used  for  the  event  distances  and  depths  encoun 
tered  in  this  research.  The  constraints  imposed  require  that  the 
amplitude  of  the  ccpstrum  peak  at  x must  be  at  least  .7  times 
the  amplitude  of  the  ccpstrum  peaks  at  both  «T  and  (a-l)x 
This  eliminates  detection  of  ccpstral  patterns  for  cases  in 
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which  the  strength  of  the  pP  arrival  is  considerably  less  than 
that  of  the  sP  arrival.  Thus,  a cepstrum  peak  at  t resulting 
from  a pP  arrival  for  an  event  will  not  contribute  to  the  CMF 
output  for  an  apparent  event  having  an  expected  sP-P  delay  time 
of  t,  without  the  significant  presence  of  a pP  arrival  for  this 
apparent  event. 

As  an  example  of  the  application  of  this  technique,  consider 
Figure  2.  At  the  top  of  the  figure  is  the  cepstrum  calculated 
from  an  event  of  known  depth  having  both  the  pP  and  sP  depth 
phases  clearly  identifiable.  This  cepstrum  pattern  is  dominated 
by  three  peaks  corresponding  to  the  sP-pP,  pP-P,  and  sP-P  delay 
times.  Upon  interpreting  this  cepstrum  for  an  unknown  event  depth, 
one  sees  that  each  of  the  three  peaks  has  the  possibility  of 
corresponding  to  the  pP  phase.  The  CMF  output  is  plotted  at  the 
bottom  of  Figure  2 and  indicates  a much  stronger  empiasis  on  the 
peak  corresponding  to  the  correct  pP-P  delay  time  f c r this  event. 
This  result  reflects  the  fact  that  this  cepstrum  pattern  primarily 
consisting  of  three  peaks  is'  in  strong  agreement  with  the  pattern 
expected  for  a single  event  having  this  pP-P  delay  time.  The 
lesser  probability  that  the  relative  location  of  these  three  peaks 
was  coincidental,  and  that  one  of  the  other  two  peaks  actually 
corresponded  to  the  correct  pP-P  delay  time,  is  indicated  by  the 
presence  of  CMF  output  peaks  at  those  delay  times  having  reduced 
amplitudes . 

Another  example  of  the  usefulness  of  the  CMF  is  shown  in 
Figure  3.  At  the  top  of  this  figum  is  a cepstrum  calculated  for 
the  same  event,  but  using  data  recorded  at  a single  station. 

The  figure  contains  four  cepstrum  peaks  having  approximately  equal 
amplitudes  making  the  interpretation  of  the  source  depth  diffi- 
cult. At  the  bottom  of  the  figure  is  the  CMF  interpretation  of 
this  cepstrum  which  reveals  a strong  emphasis  on  the  peak 


2-8 


CliPSTIUJM  MATCIini) 
F1LTI-R  OUTPUT 


IpP-P 


tine 

delay 

(seconds) 


CO 


CB 

2 


ft 


o! 


CA  <4  — 

>»  w 


pp-p 

delay  rime 
(seconds) 


l:i pure  2. 


Cepstrum  and  Cepstrum  Matched  1: j Iter  out  put 
computed  from  data  recorded  at  several  stations 
for  an  event  with  well  defined  deptli  phases. 


a m p 1 i t u cl 


CE-\f;Titu:i 


a single  station  for  an  event  with 
well  defined  depth  phases. 


corresponding  to  the  correct  plJ-r  delay  time, 
showing  the  usefulness  of  the  CMF  technique  in 
mic  source  depths  will  be  shown  in  Section  3. 


Additional  examples 
determining  seis- 


? . 3 Use  of  Seismic  Travel  Times  To  Account  For  Depth 
Phase  Delay  Time  Variations  Along  The  Coda 

In  the  previous  sections,  we  have  discussed  techniques  which 
statistically  allow  for  the  variations  in  the  differential  delay 
times  along  the  coda  in  order  to  enable  more  of  the  depth  phase 
information  to  contribute  to  the  depth  phase  detection.  However, 
allowing  for  depth  phase  delay  time  variations  through  use  of  the 
Stochastic  Cepstrum  Stacking  technique  proves  to  be  useful  up  to 
stochastic  window  widths  of  about  1 second.  Events  of  interest 
to  this  work  can  have  delay  time  variations  of  several  seconds; 
for  example  (pP-P) - (pPP-PP)  = 3.9  seconds  for  a depth  of  50 
kilometers  and  an  epicentral  distance  of  30°.  The  stochastic 
stacking  technique  could  not  encompass  such  variations  without 
a severe  loss  in  detectability  and  resolution,  unless  that 
portion  of  the  coda  containing  the  PP  arrivals  from  the  analysis 
is  eliminated  or  some  means  is  utilized  to  normalize  this  data 
to  the  direct  phase. 

f 

In  Figure  4,  we  show  an  illustration  of  such  a case.  Here, 
seismograms  are  depicted  for  recoidings  of  an  event  at  two  sta- 
tions. From  the  illustration,  we  s°e  that  the  pP-P  delay  time 
is  considerably  larger  than  the  pPP-PP  delay  time,  and,  tnere- 
fore,  the  cepstrum  peaks  cjmputed  from  samples  1 and  2 will 
not  constructively  average.  The  right  side  of  Figure  4 illus- 
trates that  the  peaks  in  the  CMF  output  would  also  lie  at  differ- 
ent delay  times  and  can  not  be  constructively  averaged  in  a 
linear  sense.  A similar  situation  is  found  with  the  data  recorded 
at  station  2,  but  here  the  delay  times  are,  in  addition,  altered 
relative  to  station  1 by  the  station  moveout.  'One  would  like  to 
be  able  to  combine  the  CMF  outputs  computed  from  the  four  samples 
in  a manner  which  would  allow  these  peaks  to  constructively  rein- 
force at  the  pP-P  delay  time.  This  is  shown  at  the  bottom  of 
the  figure. 
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If  one  calculates  the  delay  time  variations  expected 
along  the  coda  for  various  event  depths  and  distances  using 
a seismic  arrival  time  program,  the  appropriate  delay  time 
corrections  could  be  applied  to  the  CMF  output  of  each  sample 
before  averaging.  The  degree  to  which  the  CMF  output  peaks 
stacked  would  then  be  a measure  of  how  well  the  various  depth 
phase  delay  times,  computed  along  the  coda  and  at  different 
stations,  all  agreed  with  a given  source  depth  assumption. 

Using  such  a procedure,  all  of  the  depth  phase  information 
contained  in  the  seismogram  would  contribute  to  the  depth  phase 
letection . 

In  order  to  develop  such  a procedure  one  first  needs  the 
appropriate  travel  time  information.  Differential  travel  times 
such  as  sP-P,  pPP-PP,  sPP-PP.  pPcP-PcP,  sPcP-PcP,  pPPP-PPP,  and 
sPPP-PPP  are  not  easily  obtained  from  existing  tabulations. 

These  are  the  phases  which  comprise  that  portion  of  the  seismic 
coda  usually  available  for  depth  phase  analysis.  To  determine 
these  travel  times,  a ray  tracing  program  based  on  the  spherically 
symmetric  isotropic  earth  velocity  model  used  to  calculate  the 
BSSA  (Bull. Seism.  Soc.  Amer.,  Vol . 58,  No.  4(extract),  August 
1968)  seismological  tables  for  P phases  was  used  . The  program 
was  checked  by  exactly  reproducing  selected  travel  times  listed 
in  the  BSSA  tables. 

Figures  5,  6,  and  7 give  the  computed  time  differences  for 
the  (pP-P)-(pPP-PP),  (pP-P)-(pPPP-PPP) , and  (pPcP- PcP) - (pP- P) 
time  delay  differences.  These  charts  give  the  expected  shift 
of  the  CMF  output  peaks  relative  to  the  pP-P  delay  time,  com- 
puted from  portions  of  the  coda  containing  either  the  PP,  PPP , 
or  PcP  arrivals.  One  notes  that  these  shifts  are  more  signifi- 
cant for  event  depths  greater  than  -v  10  km.  Through  use  of 
these  charts,  cepstra  or  CMF  outputs  computed  from  different 
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Figure  5 

Plot  of  time  difference 
(seconds)  between  pP-P  and 
• pPP-PP  travel  time  delay 
vs.  epicenter  distance  (A0) 
for  eivon  denths. 
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Figure  6 

Plot,  of  time  difference 
(seconds)  between  pP-P  and 
pPPP-PPP  travel  time  delay 
vs.  epicenter  distance  (a^) 
for  given  depths. 
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Figure  7 

Plot  of  time  difference 
(seconds)  between  pPcP-PcP 
and  pP-P  travel  time  delays 
vs.  epicenter  distance  (A°) 
for  given  depths 


portions  of  the  seismogram  can  be  constructively  averaged  to 
produce  a prominent  peak  at  the  pP-T  delay  time.  Such  a tech- 
nique should  prove  very  useful  in  improving  both  depth  phase 
detection  and  depth  phase  delay  time  resolution. 

The  travel  time  differences  for  pPP-PP,  pPcP-PcP,  and 
pPPP-PPP  are  shown  in  Figures  8 and  9.  In  Figures  10  and  11 
are  shown  plots  of  the  pP-P  delay  times  for  the  depths  5 to 
40  km  and  5 to  100  km  respectively.  In  Figure  12,  the  PP-P, 
PcP-P,  and  PPP-P  delay  times  for  a range  of  depths  and  distances 
are  given.  The  travel  time  differences  plotted  in  Figure  12 
can  be  represented  by  single  curves  for  source  depths  between 
5 and  100  km  since  the  depth  dependence  turns  out  to  be  very 
weak.  This  points  out  that  these  phases  alone  are  not  useful 
for  depth  determination.  The  arrows  in  the  figure  indicate 
the  minimum  epicenter  distance  to  which  the  PP,  PPP,  and  PcP 
phase  exists  for  given  depths. 
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Figure  8 

Plot  of  travel  tunc  diff 
encc  (seconds)  for  pPcP- 
PcP  and  pPP-PP  vs.  depth 
(k.n)  for  given  epicenter 
distances  (A0) 
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Figure  9 

Plot  of  travel  time  differ 
ence  for  pPPP-PPP  Vs.  dept 
(km)  for  given  epicenter 
distances  (A°) 
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Figure  10 

Plot  of  travel  time  differ 
cnee  (seconds)  for  pP-P 
vs.  depth  (km)  for  given 
epicenter  distances  (a°) 
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riot  of  travel  time  differ 
once  (seconds)  for  pP-P 
vs.  depth  (km)  for  given 
epicenter  distances  (A0) 
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3.0  APPLICATION  OF  THE  NEW  DEPTH  PHASE  DETECTION 
TECHNIQUES  TO  THREE  SEISMIC  EVENTS 


Three  seismic  events,  having  known  source  depths,  were 
chosen  for  analysis  using  the  new  depth  phase  detection  tech- 
niques. These  events  had  depths  of  24,  178,  and  53  kilometers 
and  represent  a range  of  signal  to  noise  ratios  such  that 
visual  identification  of  the  depth  phase  is  easy  for  one  event, 
difficult  for  another  and  impossible  for  the  third  event. 

3.1  Analysis  of  the  Illinois  Event 

The  Illinois  Event  of  11/9/68  was  analyzed  for  differential 
travel  time  information  using  data  recorded  at  six  LRSM  stations 
(KN-UT,  MN-NV,  PGZBC,  WHZYK,  FB-AK,  NP-NT)  at  a 20  pts/sec. 
sample  rate.  This  event  was  selected  because  the  pP  phase  can 
be  accurately  identified  from  the  seismograms  (see  Figure  13) 
and  results  obtained  using  'lie  new  techniques  could  be  verified. 


3.1.1  Evidence  of  Improved  Depth  Phase  Detectability  Using 
Seismic  Coda 

The  analysis  was  carried  out  by  first  computing  cepstrums 
using  14  consecutive  12.8  second  data  samples  from  each  of  the 
seismograms.  (The  same  analysis  was  carried  out  using  samples 
which  overlapped  by  25%  and  50%  and  no  significant  changes 
in  the  results  were  noted.)  Cepstrum  phasors  as  well  as  the 
cepstrum  amplitudes  were  computed  using  the  steps  described 
in  Appendix  C.  In  Figure  14  are  plotted  the  cepstrum  stacks 


3-1 


ILLINOIS  EVENT  11/9/68 


14a 


14b 


L4c 


4d 


CEPSTRUM  STACK 

OST  - 12,3  SECOND  SAMPLE 
FROM  EACH  STATION) 


CEPSTRUM  STACK 

(5  - 12.8  SECOND  SAMPLES 
FROM  EACH  STATION) 


CEPSTRUM  PHASOR  STACK 

(5  - 12.8  SECOND  SAMPLES 
FROM  EACH  STATION) 


synthetic  ;lpsipu:i  stack 

(Tj  - 7.0;  "t2  - io.oT 


igurc  14 


CEPSTRUM  STACK  FOR  ILLINOIS  EVENT 


computed  by  averaging  six  cepstra  calculated  usirg  the  first 
12.8  second  data  sample  from  each  station  (the  start  of  these 
data  samples  was  ^ 1 second  before  the  onset  of  the  P-wave) . 

There  is  a dominant  peak  at  ^ 6 . 8 second  delay  time  corresponding 
to  a pP-P  delay  time  for  a 24  km  deep  event  recorded  at  approx- 
imately 30  . In  Figures  14b  and  14c,  are  the  cepstrum  stack 
and  cepstrum  phasor  stack  each  calculated  by  averaging  30 
cepstra  computed  using  the  first  five  12.8  second  samples  from 
all  six  stations. 

In  both  of  these  cepstrum  stacks,  a clear  peak  at  ^9.8  seconds 
has  now  emerged  and  corresponds  to  the  known  sP-P  delay  time 
for  this  event.  This  is  evidence  that  improved  depth  phase 
detectability  is  obtainable  by  including  data  from  the  seismic 
coda  in  the  analysis.  Figure  14d  is  a plot  of  a cepstrum  stack 
computed  from  data  simulated  to  have  pP-P  and  sP-P  delay  times 
of  7 and  10  seconds  respectively.  There  is  good  agreement  of 
this  theoretical  cepstrum  result  with  those  computed  from  the 
data.  (Appendix  D describes  the  steps  used  to  generate  the 
simulated  seismogram.)  Figure  15  shows  the  cepstra  calculated 
for  the  individual  stations  used  to  record  the  Illinois  Event. 
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Figure  15.  RESULTING  CEPSTRUM  STACKS  FOR  INDIVIDUAL  STATIONS 
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USE  OF  SEISMIC  CODA  IN  DEPTH  PHASE  DETECTION 


Figure  16.  Ceps t rums  computed  with  and  without  the  first 
arrival  portion  of  the  seismogram  included. 
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3.1.2  Effect  of  Later  Portions  of  the  Coda  on  Detection 


From  the  results  presented  in  Figures  14  and  16,  it  is 
clear  that  depch  phase  information  contained  in  the  first  minute 
of  coda  can  be  utilized  to  impro  e depth  phase  detection.  How- 
ever, by  including  the  next  90  seconds  of  data  in  the  analysis, 
detection  is  found  to  deteriorate  because  the  differential 
delay  times  are  substantially  different  in  this  portion  of  the 
coda.  As  was  discussed  in  Section  2.3,  one  must  deterministically 
account  for  the  delay  time  variations  associated  with  the  later 
seismic  arrivals  when  these  variations  are  greater  than  ^ 1 second. 
The  results  of  using  different  portions  of  the  coda  for  the  analy- 
sis are  shown  in  Figure  17.  In  Figure  17b  is  the  cepstrum  resulting 
from  an  analysis  which  included  the  first  3 minutes  of  the  seis- 
mograms. This  cepstrum  does  not  have  clearly  defined  peaks  and  is 
essentially  the  result  of  combining  the  cepstra  calculated  from  the 
first  minute  of  data  (Figure  17a)  with  those  calculated  from  the 
next  90  seconds  of  data  (Figure  17c)  , the  latter  having  dominant 
peaks  at  4.9  and  7.8  seconds.  The  4.9  second  delay  time  is  in 
agreement  with  the  range  of  both  the  pPP-PP  and  pPPP-PPP  delay 
times  for  a 24  km  event  recorded  over  a 20°  to  40°  range  (see 
Figures  8 and  9.)  From  Figure  12  one  notes  that  both  the  PP  and 
PPP  phases  are  present  in  the  second  minute  of  the  seismogram  for 
this  event. 

Thus,  if  one  were  to  anticipate  these  delay  time  shifts  for 
a given  depth,  one  could  then  shift  and  average  these  cepstra  such 
that  they  would  reinforce  at  the  pP-P  delay  time.  This  procedure 
would  allow  the  depth  phase  information  from  more  of  the  coda  to 
contribute  to  the  detection  and  avoid  the  deterioration  of  results 
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occurring  when  the  stochastic  stacking  window  cannot  absorb 
these  large  delay  time  variations.  For  the  Illinois  event, 
analysis  must  be  restricted  to  the  first  90  seconds  of  data  if 
one  does  not  apply  a delay  time  compensation  according  to  travel 
time  information  for  the  later  phases. 

3.1.3  Cepstrum  Matched  Filter  Interpretation  of  Result 

The  results  shown  in  Figures  14  and  16  indicate  three 
dominant  cepstrum  peaks  corresponding  to  the  sP-pP,  pP-P, 
sP-P  delay  times.  In  general  it  is  difficult  to  determine  the 
source  depth  from  this  information  since  the  relationship  among 
the  peaks  is  not  immediately  apparent  by  visual  inspection.  An 
interpretation  of  these  results  can  be  substantially  improved 
through  use  of  the  CMF  technique  discussed  in  Section  2.2. 

Figure  18  shows  the  resulting  CMF  output  plotted  with  a seismo- 
gram recorded  at  KN-UT  for  the  Illinois  event.  The  CMF  output 
shows  a single  dominant  peak  corresponding  to  the  pP-P  delay  time 
for  this  event  and  is  in  agreement  with  the  identification  of  the 
pP-P  delay  time  observed  in  the  seismogram.  This  result  reflects 
the  fact  that  all  three  cepstrum  peaks  arise  from  a single  source 
depth  having  pP-P  delay  time  of  'v  6.8  seconds,  averaged  over  the 
six  stations.  The  CMF  technique  allows  improved  and  automatic 
determination  of  the  source  depth  from  the  computed  cepstra. 
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Figure  18.  Cepstrum  Phasor  Stack,  Seismogram,  and  CMF  Output 
for  Illinois  Event. 


2 Analysis  of  the  Kamchatka  Event 


The  Kamchatka  event  of  3/15/69  was  next  analyzed  for  the 
depth  phase  delay  time  using  data  recorded  at  20  pts/sec  at 
the  four  stations  PI  IVY,  RKON,  DRCO,  and  LCNM.  The  first  minute 
of  the  seismogram  for  this  event  is  plotted  in  Figure  19.  Cep- 
stra  were  computed  using  the  first  102.4  seconds  of  data  re- 
corded at  each  station  (additional  coda  data  was  not  available 
for  this  event).  A stochastic  cepstrum  stack  was  then  computed 
using  stochastic  windows  ranging  from  0 to  3 seconds.  At  the 
top  of  Figure  20  is  plotted  the  stochastic  cepstrum  stack  com- 
puted using  a stochastic  window  of  2 seconds  which  gives  rise 
to  a significant  improvement  in  detectability  over  conventional 
averaging  techniques.  The  stochastic  window  served  the  purpose 
of  allowing  for  delay  time  variations  primarily  associated  with 
station  moveout  which  was  as  large  as  1.5  seconds  for  these  re- 
cordings. In  the  analysis  of  the  Illinois  event  the  station 
moveout  gave  rise  to  a maximum  delay  time  difference  of  * .6 
seconds  and  stochastic  stacking  gave  an  insignificant  improve- 
ment in  detectability  over  normal  averaging.  At  the  top  of 
Figure  20,  dominant  peaks  are  observed  at  * 14  seconds  and  at 
* 42  seconds-  T°  interpret  this  cepstrum  pattern,  the  CMF  tech- 
nique was  applied  and  the  resulting  output  is  plotted  at  the 
bottom  of  this  figure.  One  again  observes  additional  emphasis 
on  the  cepstral  peak  appearing  at  the  correct  pP-P  time  delay  for 
this  event  which  indicates  a source  depth  of  ^ 178  km. 
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3.3  Analysis  of  the  China  Event 


The  third  seismic  event  analyzed  for  depth  phase  informa- 
tion using  the  new  techniques  was  the  China  Event  of  6/5/64. 

This  event  was  analyzed  using  data  recorded  at  20  pts/sec  at  the 
three  stations  GG-GR,  BMO,  and  UGO.  As  is  seen  in  Figure  21, 
the  recordings  have  very  poor  signal  to  noise  ratios,  and  it  is 
impossible  for  an  analyst  to  determine  the  depth  phase  delay 
times  visually.  For  this  event,  9 cepstra  were  computed  from 
the  first  three  consecutive  25.6  second  samples  for  each  of  the 
three  stations.  These  cepstra  were  stochastically  averaged  using 
a stochastic  window  width  of  ^ 1 second,  and  the  stochastic  cep- 
strum  stack  was  computed  using  both  the  cepstrum  and  cepstrum 
phasor  techniques.  The  cepstrum  phasor  technique  gave  somewhat 
better  results  in  this  case  and  is  plotted  at  the  top  of  Figure 
22.  Two  cepstrum  peaks  stand  out  at  ^ 7 and  x 14  seconds  but 
are  not  very  convincing.  However,  by  applying  the  CMF  to  this 
data,  a more  distinct  peak  appears  in  the  CMF  output  at  the 
correct  pP-P  delay  time  for  this  event  (Figure  22,  bottom).  This 
peak  appears  at  % 14  seconds  and  corresponds  to  an  event  depth 
of  -v/  53  km  for  these  source  to  station  distances.  This  result 
gives  a striking  example  of  the  ability  of  these  new  techniques 
to  detect  the  depth  phase  and  obtain  accurate  source  depth  esti- 
mates from  recordings  of  events  in  which  conventional  analysis 
would  fail. 


Figure  21.  CHINA  EVENT  (6/5/64) 


Cepstrur.  Phasor 


Figure  22.  CHINA  EVENT  (6/5/64) 
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4.0  CONCLUSION 

During  this  research  effort,  new  seismic  depth  detection 

techniques  have  been  investigated  to  determine  whether  they 

could  significantly  enhance  one's  ability  to  obtain  accurate 

source  depths  at  teleseismic  distances.  Three  seismic  events, 

each  having  well  established  depths,  were  analyzed  using  these 

techniques  and  clear  detections  of  the  correct  pP-P  delay  times 

were  obtained  in  each  case.  For  at  least  one  of  these  events, 

conventional  depth  phase  analysis  was  not  capable  of  obtaining 

a source  depth  estimate  using  the  data  analyzed  in  this  work. 

Results  of  the  research  demonstrated  that  these  new  techniques 

represent  a major  improvement  over  conventional  depth  phase 

detection  procedures  through  advances  in  both  the  detection  and 

interpretation  of  seismic  depth  phase  information.  The  new 

analysis  method  involved  techniques  which  utilize  more  of  the 

seismic  depth  phase  information  contained  in  the  seismogram 

/ 

and  techniques  which  also  interpret  the  typically  complicated 
cepstrum  patterns  to  extract  source  depth  estimates. 

These  techniques  could  next  be  incorporated  in  a seismic 
depth  phase  analysis  package  which  would  be  of  valuable  assist- 
ance to  the  analyst.  Such  a package  would  automatically  analyze 
the  seismic  data  to  extract  accurate  source  depth  estimates  by 
determining  the  degree  to  which  the  cepstrum  patterns,  computed 
from  different  portions  of  the  seismograms  recorded  over  a suite 
of  stations,  agreed  with  those  expected  for  a given  event  depth. 
This  analysis  would  automatically  make  constructive  use  of  depth 
phase  information  contained  throughout  the  seismograms  and  should 
give  the  analyst  the  ability  to  determine  accurate  source  depth 
estimates  for  a signif icantl>  larger  percentage  of  events  then 
is  presently  possible. 
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APPENDIX  A 

STOCHASTIC  CEPSTRUM  STACKING 


A stochastic  stacking  technique  can  be  used  to  increase 
the  detectability  of  a peak  whose  position  can  vary  unpredict- 
ably  within  a certain  limited  range.  Consider  a set  of  cep- 
strums  calculated  from  consecutive  time  segments  of  a seismo- 
gram. If  the  depth  phase  delay  time  were  the  same  in  each 
segment,  then  by  stacking  (adding)  these  cepstrums,  the  ampli- 
tudes of  the  stationary  peaks  would  constructively  add,  whereas 
the  non-stationary  peaks  arising  from  origins  other  than  the 
depth  phase  delay  times  would  average  to  some  lesser  amplitude 
level.  This  is  also  true  for  the  single  cepstrum  estimate  using 
the  entire  seismogram. 

By  redefining  the  N cepstrum  values  i = 1,N  to  be 

Yi  = MAX  (x.,  j-i  - A/2,  j-i  + A/2)  , 


for  each  point  i,  and  then  by  adding  the  y arrays  one  can  in- 
crease the  detectability  for  cases  in  which  the  peak  of  interest 
moves  unpredictably  within  a window  A.  Figure  A-l  illustrates 
this  technique  using  synthetic  data.  Each  of  the  top  six  arrays 
are  constructed  of  random  numbers  having  amplitude  values  be- 
tween 0 and  1.  A peak  (indicated  by  an  arrow  in  the  figure), 
defined  by  three  points  having  amplitude  values  .8,  1.2  and  .8, 
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Figure  A-l.  Demonstration  of  Stocastic  Stacking  using 
Synthetic  Data 


is  added  to  each  of  the  six  arrays,  but  shifted  an  additional 
point  to  the  right  in  each  consecutive  array.  The  result  of 
adding  these  arrays  (Straight  Stack)  is  shown,  and  it  is  seen 
that  a clear  detection  of  the  peak  is  not  achieved  by  ordinary 
stacking.  However>  if  before  adding,  one  redefines  these  arrays 
as  has  been  outlined,  the  results  (Stochastic  Stack)  show  that 
a dramatic  improvement  in  peak  detection  is  achieved  through 
the  use  of  this  stochastic  stacking  procedure. 

Thus,  the  stochastic  stacking  technique  can  allow  for  the 
random  or  unpredictable  variations  in  the  seismic  delay  times 
which  are  not  accounted  for  when  using  travel  time  information. 
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APPENDIX  B 


STOCHASTIC  CEPSTRUM  PHASOR  STACK 


Impressive  improvements  in  seismic  depth  phase  detection, 
in  addition  to  those  gained  by  stochastic  stacking,  have  been 
obtained  by  using  information  concerning  the  phase  difference 
between  the  direct  and  surface  reflected  seismic  arrivals.  To 
see  how  this  can  be  achieved  consider  the  steps  involved  in 
calculating  a cep^crum. 


Assume  the  received  signal  F(t)  to  be  the  sum  of  the 
direct  wave  f(t)  and  a single  reflection  of  relative  amplitude 
a,  "phase"  0,  and  delay  time  t.  F(t)  can  be  written 


F(t)  = f(t)  + a [fft-r)  cos  0 + f^Ct-T)  sin  0J. 


Here  f^  represents  the  Hilbert  transformation  of  f which  cor- 
responds to  shifting  each  Fourier  component  of  f by  tt/2.  The 
bracketed  expression  then  represents  the  signal  f(t)  having 
each  of  its  Fourier  components  shifted  by  0.  This  phase  shift 
0 is  the  phase  difference  between  the  direct  and  reflected 
signal.  For  example,  if  the  reflected  waveform  differed  from 
the  direct  signal  by  only  a change  in  sign,  P would  equal  tt. 


The  power  spectrum  of  F(t)  is 
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P(w)  = p(w)  [1  + + 2a  cos  (wx-6)], 

for  w 0 and  with  p(w)  being  the  power  spectrum  of  f(t).  As 
can  be  seen  from  the  expression  for  P(w),  the  power  spectrum 
P (w)  is  modulated  by  a cosinusoidal  function  having  frequency 
t and  phase  0.  Taking  the  complex  spectrum  of  P(to)  one  would 
then  obtain  a peak  at  lag  \ having  phase  0.  Therefore,  one 
can  assign  a phase  0 to  each  point  of  the  cepstrum  calculated 
by  taking  the  power  spectrum  of  log  P(w). 


To  see  how  this  phase  information  can  aid  in  depth  phase 
detection,  consider  events  at  £ 40  km.  Here  the  direct  and 
surface  reflected  seismic  arrivals  travel  similiar  paths  once 
outside  the  vicinity  of  the  event.  It  is  then  reasonable  to 
assume  that  both  the  direct  ai.d  surface  reflected  arrivals 
would  undergo  similar  reflections,  refractions,  distortions, 
etc.,  once  away  from  the  source  region  and  that  the  phase 
relationship  established  in  the  vicinity  of  the  source  should 
be  unchanged  at  the  receiver.  If  one  considers  those  portions 
of  the  seismogram  comprised  entirely  from  P,  pP  and  sP  arrivals, 
then  the  phase  angle  associated  with  the  cepstral  peaks  re- 
sulting from  the  presence  of  the  pP  and  sP  phases,  would  be 
similar  in  all  such  portions  Other  cepstral  peaks  resulting 
from  noise,  etc.,  would  tend  to  have  more  random  phase  angle 
over  these  portions  of  the  coda.  Thus  by  treating  each  cepstrum 
point  as  a phasor  and  adding  the  cepstrums  vectorially  we  can 
enhance  the  detection  of  the  depth  phases.  Later  portions  of 
the  coda,  resulting  from  other  seismic  phases  (PcP,  PP,  etc.) 
can  also  be  included  in  this  cepstrum  phasor  stack  if  the  phase 


difference  between  P and  pP  is  the  same  as  the  phase  difference 
between  the  later  phases  and  their  associated  surface  re- 
flected arrivals  with  tt  degrees. 

These  assumptions  are  really  just  generalizations  of 
ideas  normally  assumed  about  the  direct  and  reflected  waves. 

It  is  common  to  consider  the  surface  reflected  arrivals  xc  have 
a waveform  similar  to  the  direct  arrivals  but  having  the  same 
or  different  sign  depending  on  the  angle  of  reflection  and  on 
the  signs  of  the  parts  of  the  radiation  pattern  that  radiate 
the  direct  and  the  reflected  waves.  This  is  equivalent  to 
saying  there  is  a 0 or  tt  phase  difference  between  the  direct 
and  reflected  waves.  What  we  have  done  is  to  generalize  this 
idea  to  include  all  phase  shifts,  possibly  arising  from  reflec- 
tion, refraction  or  complex  source  radiation  pattern)  from  0 
to  2tt  rather  than  just  0 and  tt  and  to  use  the  consistence  of 
0 ± tt  throughout  the  seismogram  to  aid  the  depth  phase  detection. 
The  usefulness  of  this  technique  has  already  been  demonstrated 
in  our  work  thus  far. 


APPENDIX  C 

CALCULATION  OF  THE  STOCHASTIC  CEPSTRUM  STACK 
AND  STOCHASTIC  CEPSTRUM  PHASOR  STACK 


The  following  steps  are  used  to  calculate  the  Stochastic 
Ccpstrum  Stack  and  Stochastic  Cepstrum  Phasor  Stack  for  a seis 
mogram  consisting  of  NS  consecutive  portions  each  having 
N amplitude  values: 

• Select  N sampled  (sample  rate  = 20  pts/sec) 
amplitude  values  from  the  M**1  portion  of  a 
seismogram 

XM(n) , n=0 , 1 , . . . , N-l 
and  add  N zeroes  to  interpolate  the  spectrum 

• Calculate  the  amplitude  spectrum  for  positive 
frequencies  using 


2N-1 


■ 


KtI  Vn)*-2"in:i/2N 


; j-o,i, 


N - 1 


n = 0 


with  the  Nyquist  frequency  = 10  Hz  and  i = (-1)^^. 

• Retain  only  the  lower  quarter  frequencies  of 
the  amplitude  spectrum  since  there  is  very 
little  energy  at  frequencies  above  2.5  Hz 
for  natural  events.  This  leaves  the  array 
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(AM(j),  j = 0,  1,  N/4-1) 


Remove  the  mean  and  apply  a cosinusoidal  taper 
to  the  first  101  and  last  20°6  of  the  AM  array 
giving  the  modified  array 

* 

Cam  C j ) , j = o,  1,  . , N/4-1) 

The  20 % taper  on  the  higher  frequencies  was  used 
to  de-emphasize  the  higher  frequencies 


• Add  N/4  zeroes  to  interpolate  the  cepstrum  giving 
the  array 

(Am  (j)  » j = 0,  1,  . . . , N/2-1) 

r 

At  this  stage  one  would  take  the  log  of  this 
Am  array  to  obtain  the  log  cepstrum;  but  by  not 
using  the  log  better  results  were  obtained. 


Calculate  the  Fourier  Transform  of  the  A."  array 


using 


FM^  = N7T 


N/2-1 

71  I \l  U)e'2Uijk/(N/2) 

j^O 


where  (F^(k) , k = 0,  1,  N/4-1)  arc 

complex  numbers  representing  the  positive  fre- 
quency spectrum  of  A^.  One  can  now  obtain  the 
amplitude  and  phase  of  each  ccpstrum  point  k. 

The  array 

( |F^(k)i,  k = 0,  1,  N/4-1)  is  what  we 

is  referred  to  as  the  ccpstrum  amplitude  and  the 
complex  numbers  F^(k)  are  referred  to  as  ccpstrum 
phasors . 

• To  calculate  the  stochastic  cepstrum  stack,  one 
then  calculates 

|FM(k)|  = Max  ( |FM(j)|,  j = k-A/2,  k+A/2) 

for  each  section  M and  sums  these  over  the  number 
of  sections  (NS)  used  from  a given  sci smog ram 
(a  is  the  stochastic  window  width)  . This  gives 


NS 

coo  - 2 Ifm(1>I‘wm  . 

M=1 

where  C(k)  is  the  stochastic  ccpstrum  stack  and 
is  a weighting  factor,  chosen  such  that  the 
mean  amplitude  of  each  |F^(k)|  array  is  equal. 
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• To  calculate  the  stochastic  cepstrum  phasor  stack, 
each  FM(k)  is  replaced  by  PM(n)  where  n is  the  in- 
dex of  the  Max  value  of  |F^(j)|  in  the  interval 

j-k-A/2,  k+A/2.  The  stochastic  cepstrum  phasor 
stack  is  then  defined  by 


CP (k) = 


M 


(k) 


r 
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APPENDIX  D 


SYNTHETIC  SEISMOGRAMS 


Synthetic  seismograms  are  generated  by  passing  white 
noise  through  a recursive  digital  filter  having  a hand  pass 
typical  for  seismic  arrivals  and  adding  this  signal  to  itself 
at  delays  and  i2  (corresponding  to  pP-P  and  sP-p  time 
delays).  The  recursive  digital  filter  was  designed  from 
the  following  z-transform  of  a resonator  with  poles  at 
z = r cxp(±jw  T)  and  a zero  at  q: 


H(z) 


1 -qz 


-1 


(1  - 2r  (coswrT) z' * + rZz*2) 


For  q 1 this  gives  zero  gain  at  u=0  and  a resonant  response 
with  determining  the  resbnant  frequency  and  r relating  to 
the  Q of  the  response.  We  used  q = l,  r=.9  and  ujr  = 2tt  for  the 
generation  of  the  synthetic  data  with  T being  the  inverse  of 
the  sample  rate. 
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